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Abstract 

In this paper, we present a new formulation of smoothed particle hydrodynamics (SPH), 
which, unlike the standard SPH (SSPH), is well-behaved at the contact discontinuity. The 
SSPH scheme cannot handle discontinuities in density (e.g. the contact discontinuity and 
the free surface), because it requires that the density of fluid is positive and continuous 
everywhere. Thus there is inconsistency in the formulation of the SSPH scheme at dis¬ 
continuities of the fluid density. To solve this problem, we introduce a new quantity as¬ 
sociated with particles and “density” of that quantity. This “density” evolves through the 
usual continuity equation with an additional artificial diffusion term, in order to guarantee 
the continuity of “density”. We use this “density” or pseudo density, instead of the mass 
density, to formulate our SPH scheme. We call our new method as SPH with smoothed 
pseudo-density (SPSPH). We show that our new scheme is physically consistent and can 
handle discontinuities quite well. 

Key words: hydrodynamicsmethods: numerical 


1. Introduction 

Smoothed particle hydrodynamics (SPH) is one of the methods to solve the equations of fluid 
by expressing fluid as a collection of fluid particles. It was proposed by Lucy (1977) and Gingold & 
Monaghan (1977). It is suitable for systems with large voids or systems which exhibit large structural 
changes. Thus it has been widely used for simulations of planetary science and astrophysics. 

Recently, however, it has been reported that standard SPH (SSPH) suppresses the Kelvin- 
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Helmholtz instability (Here after KHI) (Agertz et al. 2007, Okamoto et al. 2003). The reason for this 
suppression is that SSPH requires that density is positive and continuous even at the discontinuities 
of density (e.g. the contact discontinuity). This inconsistency causes large errors in the pressure 
estimate, which then causes unphysical repulsive force between particles in the low- and high-density 
region. 

So far, a number of solutions have been proposed for this problem. Examples are the use of 
pressure (energy density) instead of the mass density for the SPH formulation (we call this formulation 
density independent formulation of SPH, DISPH) (Ritchie & Thomas 2001, Read et al. 2010, Saitoh 
& Makino 2013, Hopkins 2013, Hosono et al. 2013, Rosswog 2014). Price (2008) introduced an 
artificial thermal conductivity (AC) which the pressure distribution smooths. Cha et al. (2010) and 
Murante et al. (2011) employed Godunov SPH (GSPH) proposed by Inutsuka (2002). Each of these 
solutions still has its intrinsic disadvantages. For instance, DISPH has difficulties when the pressure 
is close to zero, as is the case at the free surface, while the AC term does not exist in the original Euler 
equation. We will revisit this point in section 2.2. 

In this paper, we propose a new formulation of SPH which can in principle handle systems 
with any discontinuity. We introduce a new quantity y which we call pseudo-density, and require the 
continuity and positivity of y instead of that of the density. The quantity y follows the continuity 
equation like that of the density, but with an additional diffusion term. We can set arbitrarily value as 
the initial condition for y as long as it is positive. Except for the initial moment, y is always continuous 
and positive, as the result of diffusion, even if we give discontinuous initial distribution. Note that the 
introduction of y, which is the density of something, means we introduced effectively an extensive 
quantity associated to particle. We use the symbol Z for this quantity and we call it pseudo-mass. 
This y-Z pair is used only to construct the SPH approximation. Though they diffuse following the 
diffusion equation, this diffusion does not introduce numerical error beyond the discretization error of 
the SPH scheme, even at discontinuities. We call this formulation of SPH, smoothed pseudo-density 
SPH (SPSPH). We will show that for contact discontinuities in general SPSPH gives the results better 
than or at least as good as that of SSPH. If we choose the initial pseudo-density adequately, SPSPH 
gave the result better than that of SSPH for all tests we tried so far. Moreover it might be extended to 
handle free surface. 

The structure of this paper is as follows. In section 2, we describe the problem of SSPH and 
solutions proposed so far. We also discuss problems that have not been solved in previous studies. We 
propose our new method, SPSPH, in section 3. In section 4 we present the comparison of the results 
of test calculations with SPSPH and SSPH. Finally we present discussion in section 5 and summary 
in section 6. 
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2. Standard SPH and Its problem 

2. 1. Formulation for the standard SPH 


In SPH, the fluid is expressed by a collection of fluid particles. Physical quantities of fluid 
particles are approximated by the convolution of the quantity with the kernel function. The kernel 
convolution of a physical quantity / at the position r is defined 

/ OO 

f{r')W(\r'-r\,h)dr'. (1) 

-OO 

The function W{\r' — r\,h) is the kernel function and h is the smoothing length. The kernel 
function must satisfy the following four properties: (I) it must converge to the 5 function in the limit 
of h —y 0, (A) its integration is normalized to unity, (M) it is a function with a compact support, (IV) it 
can be differentiated at least once [C l (r) class]. 

The first order derivative, Vj-/(r), is given by 

/ OO pOO 

V r 'f(r')W(\r'-r\,h)dr'= / f(r')V r W(\r' - r\, h)dr', (2) 

by applying the partial integral and using the property (HI) of the kernel function. In order to evaluate 
the value of a physical quantity f(r) at the position of a fluid particle, we discretize Eq.(l) as follows, 


(fa)(r a ) = J2fbW(\r ab \,h a )AV b , (3) 

b 

where the subscripts a and b denote particle indices, f b is the value of / at the position of particle b 
and r ab = r a — r b . In SSPH, AV b is replaced by m b /p b , where m b and p b are the mass and the density 
of particle b. Thus f a is given by 


fa = Efb — Wab(h a ), 
b Pb 

where W ab (h a ) = W(\r a - r b \,h a ). From Eq.(2), V a / a is given by 


( 4 ) 


VJ a = Y,fb—V a W ab (h a ), (5) 

b Pb 

where the operator V Q denotes the first order differentiation by r a . In this paper, we determine h a by 
the following equation: 


h 


a 


v 



1 ID 


( 6 ) 


Here, n a is the number density of particles at the position of particle a, and D is the number of 
dimensions. The number density n a is defined as 


n a = J2 W ab(ha)- (7) 

b 

The coefficient p is a positive constant and equal to 1.6 in this paper. 

Now we derive the expression of fundamental equations of fluid for the standard SPH (SSPH). 
In SSPH the density p is obtained by substituting p for / in Eq.(4), 


3 


Pa ^ ] 7T7ij ) I l V a fc(/t a ). 


( 8 ) 


First we derive the equation of motion from the SPH Lagrangian (Springel & Hemquist 2002; 
Rosswog 2009; Springel 2010; Hopkins 2013). The Lagrangian is given by 


where 


L (q) = J2l( m b v b)-J2 m bUb + J2 x b(i)b, 

b z b b 


q = {ri,r 2 , - ■ ■ ,r b , - ■ • ,hi,h 2 , - ■ ■ ,h b , - ■ ■), 

A b 'y —1 

Ub= ^\ Pi • 


(9) 

( 10 ) 

( 11 ) 


The quantity A depends only on the entropy and is defined by A = p 1 /P, and (j> is the constraint for 
the smoothing length and equal to 

1 \ 1/D 


\n b J 


By solving the Euler-Lagrange equation for h a , we obtain 

-l 


A a = 


P a m a dp a 
Pi dh a 


1 + 


K dn a 
Dn a dh a , 


where 


C'fl'a b 


dn a 

dh n 




and d ha is d/dh a . 

Now, we can solve the Euler-Lagrange equation for r a , which is given by 
d ( dL \ dL 

dt \dv a ) dr a 

The first term is rewritten as 
d ( dL\ 

= m a v a , 


dt \dv af 

and the second term becomes 


DL ' Ppwib 9p b », h b dn b 

= ~ 2 ^ ,2 + 2^ 


dr 


b pi dr a ' ^" u Dn b dr a 

By differentiating p b and n b with respect to r a and using Eq.(7) and Eq.(8), we have 

dp b 


dr c 
dn b 
dr n 


^ ' m e ^7 a lT)) e (/j.{ ) ) (S ba d ca )y 

C 

= E V aW bc (h b )(S ba -S ca ), 


( 12 ) 

(13) 

(14) 

(15) 


(16) 

(17) 

(18) 

(19) 

( 20 ) 















where 5 ba and 5 ca are the Kronecker delta. Using V a W a b{h a ) = -V b W ab {h a ) and W ab (h a ) = W ba (h a ), 
from Eqs.(18) through (20), we have 


|^ = -£ K - n a ) V a W ab (M - £ ^ (m a - U b ) V a W ab (h b ) 
°^ r a b Pa b Pb 

Here, 12;, is the correction term derived from the variation of smoothing length as 

-l 


= 


0p a / Dn a i <9n ( 

dh 


+ 


a \ ,L a 


dh n 


( 21 ) 


( 22 ) 


(23) 


From Eq.(16), Eq.(17) and Eq.(21), we obtain the equation of motion of the following form, 

= -£ ^ K - 12 a ) V a lE- ab (/i a ) - £ (rn a - Q b ) V a W ab (h b ). 
b Pa b Pb 

The equation of energy is given by two ways. We can derive it from the law of conservation 

of energy or the first law of thermodynamics. First we present the former. 

The change of the internal energy is the same as that of the kinetic energy with an opposite 

sign. Thus the law of conservation is given by 


dm a u a 

dt 


+ 


(dm b u b \ ( d m a \v a \ 


+ 


' d m b \v b \ 
dt 2 


= 0 . 


(24) 


\ dt I \dt 2 

\ /t/ \ / a \ / u \ /a 

From Eq.(23), the change of the kinetic energy due to the pairwise interaction between particles a and 
b is given by 


d m a \v a \ 2 \ ( d m b \v b \ 2 ' 


dt 


dt 


= ~V n 


= v a ■ m a v a ) + v b • -^( m b v b ), 


Pania (m b - n a ) V a W ab (h a ) + ^ (m a - n h ) V a W ba (h b ) 


pi 


Pi 


~ Vb • (m a - n b ) V b w ba (h b ) + ^2 ( m6 _ n a ) v 6 w a 6 (/i a ) > ), 

P a m a 


Pi 


= ~Vab ■ 


Pi 


(jYl b 12 a ) a\^ab{ha) 


Dab" 9 {jYla 12 ft) ^akFfe a (/ia) j 

Pb 


(25) 


where v ab = v a — v b . From Eq.(24) and Eq.(25) we obtain 


dm a u a \ i / dm b u b \ 

V dt ) a 


dt 


+ 


^ab ' 


Pam a 

Pi 


(jn b H a )V Q l / F a &(/i a 


+v ab -^(m a -Q b )V a W ab (h b ). 

Pb 

By using constants a and /3, this equation can be expressed in the following form 


dm a u a 

dt 


P n m n 


^ab * ^ " 


P b m b 


(26) 


P 2 a 


{jTL b 12a) ^ a\Vab{ha) T V ab • [3 „ (^a 12;,) , 




( 27 ) 
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=v ab -(l-a)?^(rn b -n a )V a W ab (h a ) 

V / a P& 

+v ab • (1 - P)^ (m a - n b ) V a W ab (h b ) 


Pb 


( 28 ) 


Then we replace a with b in Eq.(28). 


' dm a u a 
dt 


= v ab • (1 - (m a - fife) V b W ab {h b ) 


+v ab -(1-P) 


Pb 
P n m n 


pi 


{jn b fi a ) VfelT a fe(/i a ). 


(29) 


Eq.(29) must be the same as Eq.(27). Therefore we have the constraint a + /? = 1. In this paper we 
chose (a,P) = (1,0). Thus the energy equation is given by 


du a \ ' -P a m a , . 

~ / J (m^ il a ) V ab • VaWafe- 


dt 


pI 


(30) 


Next we derive the equation of energy from the first law of thermodynamics. The law is given 


by 


d(m a u a ) = —P a dV a = —2^-2<ip Q 

pi 


We assume an isentropic dynamics. The equation of energy is given by 

du a rn a P a dp a 
ma dt ~ p 2 a dt' 

By differentiating Eq.(8) with respect to t, we obtain 

^ = -2>^. W ^1F 

We differentiate partially n a with respect to t using Eq.(7) to derive the second term as 


(31) 


(32) 


dn c 

dt 


^ ab ' ^ a H^afe ( h a 


By rewriting Eq.(34), the temporal differentiation of n a is given by 

-1 


dn a 

dt 


= (J2 V ab-^aW ab (h a )\ 1 + 


h a dn a dn a 
Dn a dh a dt 

of n a 

h a dn a ' 


Dn„, dh„ 


From Eqs.(22),(34) and (35) we obtain, 

dp a 


dt 


- 'y y j^lb fi a)v ab • V a Wg fe ( hg ) . 


When we substitute Eq.(36) into Eq.(31), we obtain the equation of energy. 

m °-Tr = Y — Y 2 - ( m b - fi a)v ab ■ V a W ab (h a ). 

dt b Pa 

Eq.(37) is identical to Eq.(30). Thus it satisfies the law of conservation of energy. 


(33) 


(34) 


(35) 


(36) 


( 37 ) 






















We use the equation of state for an ideal gas, 

Pa = (7 - 1 )p a u a - (38) 

2.2. Problems ofSSPH and Previous studies 

We can see that the density must be positive and continuous for the SSPH equations, 
Eqs.(8),(23) and (30) to be valid. Thus, SSPH smooths density even at the density jump of a con¬ 
tact discontinuity. As a result, the density in high (low) density region around the discontinuity is 
under (over) estimated. This error propagates to the evaluation of pressure through EOS. Eventually, 
we face a large pressure error around the contact discontinuity. This error works as an unphysical 
surface tension, resulting in the suppression of fluid instabilities (Agertz et al. 2007). To overcome 
this difficulty, several modifications of SSPH have been proposed. 

Price (2008) proposed the use of AC. In his approach, AC makes the energy smooth so that 
the pressure distribution can be flat. Thus both the density and the pressure become smooth. It works 
fine if the discontinuity is due to the jump in the thermal energy. However, it is not clear how we can 
handle the discontinuity in the chemical composition. 

Price (2008) and Read et al. (2010) showed that the KHI takes place with the smoothed 
pressure formulation which was developed by Ritchie & Thomas (2001). Saitoh & Makino (2013) 
clarified the mathematical and physical implication of the scheme used by Ritchie & Thomas (2001) 
by means of a volume element AV a given by m a u a /q a . The quantity q a is the energy density defined 
by Qa = p a u a , and it is proportional to the pressure in the case of the ideal gas. Thus, in DISPH, 
pressure should be positive and continuous. DISPH can deal with the contact discontinuity without 
any difficulty. However it still has some weak points. For example, it cannot evaluate the volume 
element of the fluid particle of which the pressure is very small (e.g. near the free surface of water). 

Cha et al. (2010) and Murante et al. (2011) adopted GSPH, which greatly improves the 
pressure wiggles even in the shock tube test involving strong shock with the mach number of 10 5 . 
GSPH was originally proposed by Inutsuka (2002). They showed that the KHI grows well due to 
the good behavior of GSPH at contact discontinuities. It is, however, difficult to handle a non-ideal 
gas with GSPH, since one needs to solve a non-ideal Riemann problem and it is computationally 
expensive. 

Garcia-Senz et al. (2012) introduced an integral approach in order to evaluate the first deriva¬ 
tive of the SPH approximation. This method can reduce the error of the gradient. With this method, 
they showed that the KHI grows even though they used the SSPH. Recently, Rosswog (2014) com¬ 
bined this method with DISPH. Although this method is very efficient, it is unclear how to deal with 
problems involving free surface. 

Ott & Schnetter (2003) proposed yet another method. In their method the volume element 
AV a is defined as 1 /n a , where n a is the number of density of particles. It is defined as 

n a = Y,W ab {h a ). (39) 

b 
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This method performs well when the density evolution is continuous. If this is not the case, for 
example if two stars merge, this method would perform poorly, because the distribution of n a contains 
discontinuity. 

We extend the idea of DISPH so that we can handle the region where pressure is quite low. 
Here, we introduce a new quantity which is a virtual quantity and is not a physical one in order to 
evaluate a volume element. In the following, we describe this new formulation of SPH. 

3. New SPH 

In this section, we present the formulation of our new method, SPSPH. In SPSPH, we use the 
pseudo density y and its associated pseudo mass Z, to obtain the volume element AV. In addition, we 
let y diffuse, following the diffusion equation to guarantee that its distribution is (or will be) smooth 
everywhere. 

3.1. Equation for the pseudo density 

We derive the time evolution equation of the pseudo density y with diffusion. First we drive 
the diffusion equation in the Eulerian view. We define the quantity j as the flux density of Z by 
diffusion, 

We derive j which satisfies the following property. The quantity j depends on Vy, because it 
should reduce the jump in the distribution of y. The simplest form is 

j = -AiifVy. (41) 

The coefficient Aif should be positive in Eq.(41). The diffusion equation of y is 

= D dii V 2 y + VAuf • Vy. (42) 

dif 

In this paper we use Aif which is spatially constant. Thus the second term vanishes. In section 3.6, 
we discuss Aif which depends on r or y. 

Since y should satisfy the continuity equation, the time evolution equation should contain the 
advection term as 

^ = AifV^-V-M. (43) 

We then derive the equation for y in the Lagrangian view. The Lagrangian derivative d/dt is 
given by d/dt = d/dt + v ■ V. Therefore the equation is given by 

= D dii 'V 2 y -yV-v. (44) 

dt 

The first term in the right-hand side of this equation indicates the evolution through diffusion and the 
second term indicates the change of pseudo density through compression or expansion. 
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3.2. Generalized volume element and SPH equations 


In this section, we derive the set of fundamental equations of SPH based on the pseudo density 
y. We start with the volume element A V a given by 

AK = (45) 

Va 

where the pseudo density y a is given by 


y a = J2 Z bWab(h a ). (46) 

b 

Here, the quantity Z is an extensive quantity associated with y. In other words, y is the spatial density 
of Z. We call Z psuedo mass and y pseudo density. 

First we derive the equation of motion from the SPH Lagrangian. The Lagrangian is given by 


L (<7) = Z o (' m b v b ) - Z m b U b + Z X b(j)b, 

b Z b b 


(47) 


where 


q = (r 1 ,r 2 ,---,r b ---,hi,h 2 ,---,h b ,---), 
_ A b (m b y b \i~ x 
Ub ~ 7 - 1 V Z b ) 


(48) 

(49) 


The function </> is the constraint for the smoothing length and is the same as Eq.(12). 

As in section 2.1, by solving the Euler-Lagrange equation for h a and r a , we obtain the equation 
of motion in the following form 


m. 


iV a — Z 


Pa Za 

vl 


(z b -n a )v a w ab (h a 


Y,^(z a -n b )v a w ba (h b ). 

b 


Vb 


Here, Hf, is the correction term derived from the variation of smoothing length as 

-l 


o,= 


dy a ( Prig dug 
dhg \ h a dhg 


(50) 


(51) 


where 

^ = Z Zbd K W ab (h a ). (52) 

dhg Y 

We obtain the equation of energy from the equation of motion again as in section 2.1. It is 
given by 

rhi P 7 

-37 = Z ( Z b ~ ^a) v ab • V a W ab {h a ). (53) 

dt V Va 

In the case of an ideal gas, pressure P a is given by 


Pa 


(7-1) 


^^a^aVa 


(54) 


This formulation satisfies the law of conservation of energy. In addition, it conserves the linear 
momentum, the angular momentum and the total mass. The temporal differentiation of the total linear 
momentum is given by 
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d_ 

dt 


E m ^ = E 


dv n 


m n 


dt 


= -EE 


PaZ a 

vl 


(. Z b - n a ) v a w ab {K ) + ^( Za - n b ) v a w ab (h b 

Vb 


( 55 ) 


The sign of the right side changes, while the absolute value does not change, when we exchange the 
indices a and b, because the indices a and b are anticommutative in V a W 0 &- Thus the value of the 
right side is zero and it shows that our formulation conserves the linear momentum. The temporal 
differentiation of the angular momentum is given by 


d x , 

— 2^r a x m a v a 


= E r « x 

a 

= -EE 

a b 


m, 


dv a 
1 dt ’ 
PgZ, 
Va 


(yZ b T abPab(ha 


b 2 b (Z a ~ Qb) r a bFab(h b ) 
Vb 


(56) 


where the function F ab satisfies r ab F ab = V a W ab and r ab equals r a x r b . The indices a and b are 
commutative in F ab and anticommutative in r ab . Thus the value of the right side is zero and it shows 
that our formulation conserves the angular momentum. In SPSPH we express the density of the 
particle a as p a = m a y a /Z a . The total mass is given by 

WlbVb Z b 




b ^ b 

= E m &- 


Vb 


(57) 


Therefore, SPSPH conserves the total mass because particle mass is the constant with respect to time. 
3.3. Implementation of the diffusion term 

First, we derive the SPH expression of the first term in the right side of Eq.(44) for SPH. The 
SPH expression for the Laplacian has the following form (e.g. Brookshaw 1985): 

t'ab ' ^ afFabiflf) 


V 2 A a = 2Y / (A a -A b )AV b - 

b 

Therefore, the diffusion term of y is given by 


| fab | 



= 2 D dii J2{ya-Vb) 

dif b 


Z b T ab ‘ V (A^abiftd) 


yb 


\ r ab\ 


(58) 


(59) 


What we actually need is the equation for the pseudo mass Z, not for y, since we obtain y from 
Z using Eq.(46). In the following, we derive the diffusion equation of y for Z. The requirement is 
that diffusion doesn’t change the Lagrangian. Therefore Z must evolve so that the volume of particle 
is not changed by the diffusion of y, because Lagrangian L{y , Z) has the form L(Z/y). Therefore, 

lZ \ i / a/ \ / / a,v \ 

= 0. (60) 


' d 3 \ 


dt 


dif 


y \dt J 


dif 


y 
















Thus the diffusion equation for Z is expressed as 

f dZ a ^ _Za(ctyo\ Z, 

dt 


dt 


= 2D di {—J2(y a ~y b ) 

dif d a b d b 


Z b fab ‘ V aW^abiha 


- dif y> 

and it gives to the equation of time evolution of Z. 


r ab 


dZ, 


= 2D di{ -J2(y a -y b ) 
dt y a b Vb 


Z b 'Tab ' V a'Wabi.h'a) 


\ r ab\ 


( 61 ) 


(62) 


3.4. Time step criterion and stability analysis 

The time step for integration is limited by the Courant condition for numerical stability. Since 
we use the leap-frog method for time integration, the time step A^cfl is given by 

At C FL = min dt a , (63) 

a 

2 tin 


dta t \ 


CFL 


maXbOf 


(64) 


It is used as the time step for SPH without diffusion. In this paper, we set C C fl = 0.6. We need to 
derive the maximum time step Atdif for the diffusion equation of SPH, because it can become smaller 
than At C FL- 

We derive the time step for the diffusion equation of SPH by the linear stability analysis of 
Eq.(61). We perturb the pseudo density y and the pseudo mass Z from the uniform state which 
satisfies y a = y° : Z a = Z° and h a = h° for all particles a. The perturbations of y and Z are defined as 
5y and 5Z. By perturbing Eq.(61), it is given by 


d(Z a + 5Z a 
dt 


0 n Z a + 5Z a ^ .Z b + SZ b r ab -V a Wab(h a 

= 2-Ddif-—=— 2^{y a — Vb + oy a - dy b ) -—-;--- 

y a + oy a V Vb + oyb 


r ab 


(65) 


We start analysis by y a = y°,Z a = Z° and h a = h° for all a. Thus 


d(Z° + 5Z a ) 
dt 


on Z° + 5Z a Z° + SZ b r ab -V a W ab (h °) 

ZDa ^W»V y ~ v +Sv °~ Svb) ZT^ —FS— 


The perturbation equation to the first order is expressed as 

d(5Z a ) on .Z 0 r ab -V a W ab 

-= ^D dd —}^(dy a -dy b ) — - 


dt 


fab 


( 66 ) 


(67) 


r b v 

We can express W ab (h°) as W ab for convenience, because h° is the same for all particles. In the 
following, we consider the simplest case of a ID problem in which particles are placed in equal 
spacing. In this case, 5Z a can be expressed by the fourier series as follows 


<5Z a (r) = £Ae"V^ 


( 68 ) 


where a; is a complex number. Consider the perturbation of Z with the wave number k. 
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(69) 


5Z a , k (r) = Ae ut e ikra . 

The perturbation of y is then given by 

5y a , k (r) = J2A m e^e^W ac , 

C 

and Eq.(67) with the wave number k becomes 

Z° 


f e u(t+At di t) _ e ut\ e ikr a 


At, 


dif 


= 2B *>jEE 


e ut e ikr c 


W n 


J2e ut e ikrc W t 


(70) 

r ab -VW ab 


be 


if 


ab 


(71) 


We divide the equation by exp(cut) exp(ikr a ) and use that the particle spacing equals Z°/y°. Thus, 
we obtain 




At, 


dif 


Z° 

1° 

z° 


= 2A U(3 £ £e 


ikr c 


y° b 


'W ac -J2e ik ^»-^w bc 


z° r ab -vw c 


ab 


yU 


ab 


= 2D^Y,\H w )-Hw) e ~ ikn ' 

D —ikr ba 


r 


r ah -vw ( 


ab 


ab 


2 ZWW- ir £(l 

y h 


fab ■ VW, 


ab 


(72) 


ab 


where J-(f) expresses the fourier transform of a function /. We used W ac = W ca . Note that Eq.(72) 
converges to zero under a = b in the real space in the limit of r ab —» 0. The lowest order of the function 
V aWab respect to r ab is 0(0). Thus, that of r ab ■ VW ab /r 2 ab is 0(r~ b ). If [1 — exp (ikr ab )\ converges 
to zero faster than 0(r ab ), Eq.(72) converges to zero. From THopital’s rule, we obtain 


Re 


lim 

r a b-+ o 


1 - exp (ikr ab ) 


— sin kr ab 

= hm —— 7 - 7 — = (J. 
r a b^-0 69(0) 


(73) 


0(r\ b ) 

Thus, Eq.(72) converges to zero for r ab —> 0 in the real space. Therefore, by defining the function g as 


g(r a b) = 


E c 


VWa 


r ah -Vw ah 


(b = a) 
(b^a) 


we can obtain the following equation: 


a tdAi dif 


At, 


dif 


= 2D dii T(g)T(W). 


(74) 


(75) 


The condition that the perturbation damps is a positive At dif that satisfies | exp(ccAt d if)| < 1 
exists for all possible values of k. Therefore the stability condition is 

-1 

0 < Af di f < min ■ 


T'D m F(g)F(Wy (76) 

The upper limit of the wave number /;: rn . )X is 2ir/ Ax where Ax is the particle spacing. From 
Eq.(76), the function T(g)T(W) must be negative for all possible values of k. First, we show that 
F(g) is negative. The kernel function W(\r ab \) takes the maximam value at \r ab \ = 0, and is decreas¬ 
ing for both sides. Therefore, the gradiant is positive for r ab < 0, and negative for r ab > 0. All terms 
in the summation of g(r ab ) is negative, and therefore g(r ab ) is negative. Tig) is given by 
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( 77 ) 


Ha) = ^g(r ab )z itr “, 

y u 



r a b ■ VW ab 


r 


2 

ab 


The quantity 1 — exp (ikr ab ) is positive, whereas F(g) is negative. Thus the function F(W) must be 
positive for the possible k. In this paper we use the Wendland C 2 and C 4 function for one dimension 
and two dimensions, whose fourier transform are positive for all wave number 0 < k < /c max as kernel 
function, as proposed by Dehnen & Aly (2012). We note that the perturbation does not grow within 
a finite time even though the fourier transform of these function can become non-positive for an 
inhomogeneous particle distribution, since dZ/dt is small enough. 

The maximum value of Atdif can become very small. Therefore, we consider the option to 
use separate time steps for the diffusion equation and the hydrodynamics. We use At CF L for the 
hydrodynamics and use Af^ p for the diffusion in practice. The time step Af'ljj 1 ’ is given by 
AtcFL 


Afd7 = 


M 


where, 

M — 
We use At dif as 

Afdif = 


At 


CFL 


At 


dif 


(min a 6 |r a6 |) s 


(78) 


(79) 


( 80 ) 


(279 dif) ’ 

It is easy to show that this is the maximum time step for the equally spaced particles in the one 
dimension case. Dehnen & Aly (2012) showed T{W) < 1. From Eq.(74), 

r ab -VW ab 


m< 2 £ 


-2£ 


f/ u 

Z°w ab 


ab 


b y 
2 


o 


ab 


< 

mm bK b 

Thus the maximum value of the function T(g) respect to r ab satisfies 

2 


maxj r ((yf) < 


min abr 2 ab 


( 81 ) 


(82) 


and we have Eq.(80). 

Figure 1 shows that the timestep obtained using of Eq.(80) is smaller than that using of Eq.(76). 
Thus, the time step derived from Eq.(80) satisfies the condition of Eq.(76), independent of the value of 
the diffusion coefficient, because both equations depend on inverse proportion of diffusion coefficient. 
Therefore we use the time step derived from Eq.(80). 
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Fig. 1. Time steps derived from Eq.(76) and Eq.(80). The solid line shows time steps evaluated by Eq.(76) and the 
dashed line represents those evaluated by Eq.(80). 

3.5. Artificial Viscosity 

To deal with shocks, we need to add an artificial viscosity term to the energy equation and the 
equation of motion. The viscosity term for the equation of motion is 

[VW ab (h a ) + VW ab (h b )} 


dm a v a 

dt 


= -^m a rn h F ah Il ab - 
AV b 


and that for the energy equation is 


dm a u a \ 1 \ ^ 

-I = - 2^ m a m b F ab U ab v ab ■ 

AV Z b 


[VW ab (h a ) + VW ab (h b )} 


(83) 


(84) 


V dt 2 V "'"'’” -- 2 

Here, U ab gives the viscosity and F ab is a “switch” function to reduce shear viscosity (Balsara 
1995). There are several different forms of the function II ab . In this paper we adopt the function 
proposed by Monaghan (1997). It is expressed as 


Hah = 


-a 


Slg 

Vab^ab 


Ft{ fab ‘ fab) i 


(85) 


{pa + Pb) ' 

where H{x) is the Heaviside function, v = cs a + cs b — 3u) ab , u ab = v ab ■ r ab /r ab and cs is the sound 
speed. In this paper we set a vls = 1. The mass density p is calculated by p = m/ AV. We used the 
standard Balsara switch (Balsara 1995) for F ab . It is given by 


Fab = 2 {Fa 
F n = 



|v- v a \ 


|v-« a | 

+1 

v x v a 

| + e b cs a /h a 


( 86 ) 

(87) 

where e b is a small constant and in this paper we set e b = 10 -4 to prevent numerical overflow. The 
divergence and rotation of velocity are 


v-v a = -'£ l —v ab -vw ab (h a 
V Va 


( 88 ) 
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(89) 


7 

V X V a = -J2~ V ab X VW ab (h a ). 
b Va¬ 
in SSPH, we used density p and mass m in place of y and Z. 

3.6. Possible choice of the diffusion coefficient 


So far we assume that the diffusion constant, D d a is actually constant in space. However, it is 
possible to change D di {, since the diffusion itself does not introduce numerical error. In the limit of 
Ddif —> 0, our new SPH scheme is reduced back to SSPH. Thus, where SSPH works fine, or in other 
words, after y has become reasonably smooth, it makes sense to reduce D d[i so that we can increase 
At dif . For example, consider a polynomial form expressed as 


D dii = J2B n \Vyff 


(90) 


where B n and //,, are positive actual constants. Then Eq.(42) becomes 

(dy 

\ dt 


dif 


Y,B n \Vyf n ~ l (\Wy\W 2 y + p n V\Vy\ Wy) 


(91) 


However, in this paper, we use the diffusion coefficient D dd which is constant in space and 
variable in time. One simple choice of D d[i is such that At dif = A/ CF l- Such D d[ f is given by D M 
which is defined as 

° M = iAW <92) 

where Ax is the minimum particle spacing. In this paper, we use the initial particle separator as 
Ax, since the density change in our tests is small. If one uses this method for more realistic cases, 
it is necessary to adopt the minimum particles distribution at each time step and/or to evaluate the 
diffusion coefficient for each particle. For example, 

°m = (93) 

ZZAtcFL 

where A is a positive constant and equals to about 2.8 rj with the Wendland C 4 kernel for two dimen¬ 
tions. 


4. Numerical tests 


In this section, we compare the results of SPSPH to those of SSPH. In all tests, we found 
SPSPH gives better, or at least similar, results, compared to those of SSPH. 

In section 4.1, we show the results of the “Square test” proposed by Saitoh & Makino (2013), 
in which the evolution of a square-shaped high-density fluid embedded in a low-density fluid is solved. 
We then show the results of the shock tube tests in section 4.2 and those of the KHI tests in section 4.3. 
We give the results of square tests with an extreme density contrast in section 4.4. Finally, we show 
the results of the one-dimensional hydrostatic equilibrium tests and discuss the remaining difficulty 
of SPH induced by inhomogeneous particle distribution. 
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4.1. The square test 


This test was first used by Saitoh & Makino (2013). The computational domain is a two- 
dimensional square of unit size, —0.5 < x < 0.5 and —0.5 < z < 0.5, with the periodic boundary 
condition. The initial density is given by 

f p = 4 — 0.25 < x < 0.75 and 0.25 < £ < 0.75, 

I p = 1 otherwise. 

We set 7 = 5/3 and initially P = 2.5, v x = v z = 0. 

We used two different initial distributions for y. In the first one, the initial value of pseudo¬ 
density y init is unity everywhere. The second one is that the initial pseudo-density is the same as the 
initial density, i.e.,y lD it = p. In this way, we can see if our scheme can handle initial discontinuity of 

V- 

We expressed the density difference in two ways. In the first one, particles in the two regions 
have the same mass and the spacing between particles is changed. Thus the number density of particle 
is different in the two regions. The particle mass in this case is set to 1/16276. In the second one, 
the particle mass is changed so that the particle spacing can be the same in the whole region. The 
particle mass for the high density region is 1/4069 while that for the low density region is 1/16276. 
We summarize our runs in table 1. 

Table 1. Summary of the runs for the square test. 


Model 

SPH scheme 

Particle mass 

2/init 

D<nf 

run 1 

SSPH 

equal 

- 

- 

run 2 

SPSPH 

equal 

l 

Dm 

run 3 

SPSPH 

equal 

density 

10 D m 

run 3 a 

SPSPH 

equal 

density 

Dm 

run 4 

SSPH 

unequal 

- 

- 

run 5 

SPSPH 

unequal 

1 

Dm 

run 6 

SPSPH 

unequal 

density 

Dm 


For run 3 we let y diffuse for 0.1 unit time with the diffusion coefficient of 500(Ax) 2 , before 
we start the calculation. Here, Ax is the initial particle separation in the high density region. In 
addition, we set D di{ = 10 D M for run 3 as is seen in table 1. This is because the diffusion is insufficient 
in the case of D dli = D M (run 3a) and run 3a gives an inadequate result. One might consider that the 
use of 10 Dm results in a high computationally cost. However, the actual computational cost of run 3 
is only twice as expensive as that of run 3a. We consider that this increase of the computational cost 
is much smaller when we simulate complicate scientific phenomena in realistic situations. 

Figures 2 and 3 show the time evolution up to t = 8. It is clear that SPSPH can handle the 
contact discontinuity quite well while SSPH cannot. SSPH made the square region almost completely 
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circular. 

Figures 4 and 5 show the cross section of the pressure and the acceleration for SSPFI and 
SPSPH. The SSPFI results show large errors at boundary. These errors in SSPH work as the unphys¬ 
ical surface tension (Saitoh & Makino 2013). As a result, SSPH cannot maintain the initial particle 
distribution (figures 2 and 3). 

The results with SPSPH are far better than those with SSPH, as expected, even when y is 
initially discontinuous. The reason for this behavior is that it requires the continuity of y instead 
of that of the density. Even when y is initially discontinuous, SPSPH can make it continuous and 
smooth. Figure 6 shows the cross section of the density (SSPH) and the pseudo density (SPSPH). 
We can see that pseudo density y has become pretty smooth in SPSPH, while the jump in the density 
is well maintained. It indicates that SPSPH performs well even when the initial distribution of y 
contains discontinuity. 


4.2. Shock tube tests 
4.2.1. Sod shock tube test 

The initial condition of one-dimensional shock tube test is the same as that used in Sod (1978). 
The computational domain is —0.5 < x < 0.5 with a periodic boundary condition. We placed the 
discontinuity at the origin by setting initial condition: 

f P = 1 P = 1 v x = 0 x <0, 

1 p = 0.25 P = 0.1795 v x = 0 x > 0. 


We actually used a “smoothed” initial condition. The smoothed initial pressure is given by 


P = 1 
P = 


‘■[(C'x) 3 — 3Cx\ + 


Ph+Pi 

2 


X<-C’ 

±<x< p , 


(96) 

P = 0.1795 x >h- 

The coefficient C is an arbitrary constant and we used 103.22 in this paper. I), and /) are pressures 
in the high and low density regions (i.e., 1 and 0.1795), respectively. The smoothed initial density 
distribution is given by 


x < — 

Ph.+Pi _ j_ 


C’ 


P = 1 

p — £cz£h [(Cx) 3 — 3Cx] + enczEL < x < ^ 
p = 0.25 X > 

where p h and p t are densities in high and low-density regions (i.e., 1 and 0.25), respectively. 
The position of the particle a, x a satisfies 

rx a 

/ p(x)dx = m, 

1 X a -1 


(97) 


(98) 


where, m is the mass of a particle, and we used m = 1/1600. We carried out two cases, i.e. y m \ t = 1 
and 2/i n it = p. The initial pseudo mass is given by Z — my/p = m/p for y init = 1, and Z = m for 
Pinit = p. We set A,* = D m in this test. 
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Fig. 2. Results of the square test, for runs 1 to 3 in table 1 from top to bottom, snapshot at t = 0.1,0.5,1.0,8.0 are 
shown. In this test the density ratio is 4:1 with the mass ratio of 1:1. The top row is the result of SSPH, while the 
middle and bottom rows are results of SPSPH in which y irlit is 1 and density, respectively. 
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Fig. 3. The same as Fig. 2 but for runs 4 to 6 and with the mass ratio of 4:1. 
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Fig. 4. Pressure (left) and acceleration (right) in the x direction for three different models (rani,2 and 3). Physical 
quantities which are only in the region of \y\ < 0.05 at t = 0 (crosses) and t = 0.1 (circles). 
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Fig. 5. The same as Fig. 4 but for runs 4 to 6 and with the mass ratio of 4:1. 
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X X 

Fig. 6. The cross-section of the density and pseudo density. Values of particles in the region of \y\ < 0.05 and 
plotted. The time t = 0.1. The left hand-side panel shows the results of run 4 in table 1 (SSPH), and the right 
hand-side panel shows that of run 6 in table 1 (SPSPH). For the SPSPH run, we show both y (crosses) and density 
(dots). 

In Figure 7 we can see that the behaviors of SSPH and SPSPH calculations are very similar, 
even though the volume elements used are quite different. Thus, we can conclude that the use of 
pseudo density with artificial diffusion does not affect the behavior of SPH scheme for the standard 
shock tube test. In this case, the amplitude of the pressure wiggles of SPSPH with y init = 1 are 
comparable with or smaller than those of SSPH. In other words, SPSPH is not much better than SSPH 
in this case. The reason is that both SPSPH and SSPH have the error caused by the discontinuity of 
particle spacing. We explain this further in section 4.5. If we choose an adequate y in it , the absolute 
value of the jump in pressure can be reduced to about half of that obtained in the SSPH run. This 
is because the distribution of the fundamental quantity y in the SPSPH run with y ini) ,= p is much 
smoother than that with y init = 1 (see figure 8). 

4.2.2. Strong shock tube test 

Here we present the result of a strong shock test, similar to what is performed in Toro (2009). 
The computational domain is —1.0 < x < 1.0 with a periodic boundary condition. We placed the 
discontinuity at the origin by setting initial conditions as 

f p=l P = 1000.0 v x — 0 x <0, 

{ ' ’ ( 99 ) 

| P = 1 p = 0.01 v x = 0 x>0. 

We used ci;av = 2 and the same smoothed distribution as in Eq.(96) with C = 129.025, y init = 1 
everywhere, and m = 1/500. We set D^f = Dm, however M equals 7 with SPSPH. 

In Figure 9, again, SSPH and SPSPH behave similarly. There are some differences. For 
example, the overshoot in the density at the contact discontinuity is smaller for SSPH, while the 
jump in the pressure is smaller for SPSPH. For shock tube tests, there is no reason to expect big 
improvement over SSPH. 
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X X 

Fig. 7. Result of the sod shock test at t = 0.1. Panels in top, middle, and bottom rows show the results of SSPH, 
SPSPH with t/i n it = 1, that with j/j n j t = p, respectively. The left hand-side panels show the density and the right 
side-hand panels show the pressure. Black circles show the numerical results, and red curves analytic solutions. 
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Fig. 8. The distribution of the pseudo density y for the Sod shock tube test. The left and right panels show the 
results of SPSPH with y\ n ix = 1, and with y; n i t = p, respectively. 


Q. 


Q. 





X X 


Fig. 9. Results of the one dimensional strong shock tube tests at t = 0.012. Panels in the top and bottom sides show 
the results of SSPH and SPSPH runs, and left and right panels show the density and pressure. Black circles show 
the numerical results, and red curbs analytic solutions. 
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4.3. Kelvin-Helmholtz instability tests 


The KHI test is useful to investigate the ability of SPH schemes to handle hydrodynamical 
instabilities, (e.g., Okamoto et al. 2003, Agertz et al. 2007, Price 2008) 

We performed two-dimensional calculations and use the computational domain that is a square 
of an unit size, —0.5 <x< 0.5 and —0.5 < z < 0.5, with a periodic boundary condition. We make 
the contact discontinuity by setting initial conditions as 


p = 2 — 0.25 < z < 0.25, 
p = 1 otherwise. 


( 100 ) 


We set P = 2.5, 7 = 5/3 and v x (= v x> h) = 0.5 in the dense region ,v x (= v x> i) = —0.5 in another. The 
initial velosity perturbation is 


v z = — -—-A sin ■ 


2.07r(a: + 0.5) 


(0.225 < |*| < 0.275) 


where A = 1/6 and A = 0.025. The growth timescale of the KHI is 

Kph +pi) 


Ach = 


( 101 ) 


( 102 ) 


yjPhPl .h ^x,l\ 

For our setup, tkh — 0.35. This setup is the same as that used in Price (2008). 

We used two different initial distributions for y. In the first one, the y init is unity everywhere. 
In the second one, y init = p. In this way, we can see if our scheme can handle the initial discontinuity 
of y. For this run, we used D dif = D M and At dif = A£cfl- Unlike in the square test, we did not let 
pseudo density diffuse before we started the calculation. Particles in the two regions have the same 
mass. Thus the number density of particle is different in the two regions, 262144 in the dense square, 
and 131072 in the other region. 

Figure 10 shows the time evolution up to t = 8tkh- It is clear that SPSPH is much better than 
SSPH in dealing with KHI. With SSPH, the perturbation grow but the roll-like structure characteristic 
of the KHI is suppressed. Moreover rolls break apart by t = 4 — 8tkh• These are due to the effect of 
the artificial surface tension at the boundary of two fluids. In two SPSPH runs, the KHI grows well, 
and there is no effect of the artificial surface tension. SPSPH can handle hydrodynamical instability 
even if y is initially discontinuous. 


4.4. The square test with extreme density difference 

Here we present the results of the square test as in section 4.1, but with a much larger density 
construct. The density of the high-density region is 100 instead of four. We only did the cases with 
different mass particles since with equal-mass particles, the difference in the particle number density 
would become too large. For the SPSPH run with y init = p, we let y evolve for 0.01 time unit before 
we start time integration. Figures 11 and 12 show the results, we can see that SPSPH can deal with 
very large density contrast much difficulty, even when the pseudo density is not initially continuous. 
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P 

Fig. 10. Results of the KHI tests with the density ratio of 2:1. Density distributions at t = 1,2,4 and 8 are shown 
from left to right. Panels in the top, middle, and bottom rows show the results of SSPH, SPSPH with y; n it = 1, that 
with r/init = p , respectively. 
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Fig. 11. Results of the square tests with the density ratio of 100:1. Density distributions at t = 0.1,0.5,1.0 and 8.0 
are shown form left to right. Panels in top, middle, and bottom rows show the results of SSPH, SPSPH with t/i n it=l, 
that with y init = p, respectively. 
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Fig. 12. The same at Fig. 6 but for the density contrast 100:1. In the right-hand-side panel, squares, crosses and 
dots show the profile of y(t = 0), y(t = 0.1) and density(f = 0.1), respectively. 

4.5. The one-dimensional hydrostatic equilibrium test 

Here, we investigate the ability of SSPH, DISPH and SPSPH to handle the contact discon¬ 
tinuity which exists in a hydrostatic equilibrium state with a discontinuous particle spacing. In this 
test, we consider the one-dimensional computational domain —0.5<x<0.5 with a periodic bound¬ 
ary condition. This computational domain is filled with a gas of 7 = 1.4. The initial condition is as 
follows: 

p = 1 P = 2.5 v x = 0 x < 0, 
p = 0.25 P = 2.5 v x = 0 x > O.fhn 

The number of particles is 100, and mass of each particle is 1/160. Thus, the interparticle distance is 
1/160 for x < 0 and 1/40 for x > 0. The initial value of pseudo density y ini{ is set to unity and that of 
Z in x < 0 is 1/160 and that in x > 0 is 1/40. The diffusion coefficient for SPSPH is set to D dii = D M . 

In the left panels of figure 13, we show the distributions of pressure and acceleration for the 
one-dimensional test for SSPH, DISPH, and SPSPH (t = 0). We can see that the distributions of 
pressure and acceleration for SPSPH are identical to those for DISPH. However these for SSPH are 
different. Although the wiggles of acceleration for SSPH is comparable to these for SPSPH and 
DISPH, that of pressure for SSPH is more than three times larger than those for SPSPH and DISPH. 
Note that, of course, these quantities should be uniform throughout the computational domain. Hence, 
these wiggles are induced by the asymmetric distribution of particles. 

With these results, we can now understand the results show in section 4.1 and 4.2 better. In the one¬ 
dimensional shock tube tests, we see the weak pressure wiggles at the contact discontinuities. Thus 
wiggles are caused by the inhomogeneous particle distribution. Even in DISPH, we observed these 
weak wiggles (see figure 1 in Saitoh & Makino 2013) at the contact discontinuity. In the panels in 
the right-hand side of figure 13, we show the result of the same hydrostatic test as in the left-hand 
side panels, but in two-dimensional calculation. The density contrast is 1:4 in both cases, but the 
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ratio of interparticle separation is 1:4 in ID and 1:2 in 2D. Thus, we can expect that the effect of this 
inhomogeneity in the particle distribution would be smaller for 2D (or even so for 3D) calculations 
than in ID calculation, even when the density contrast is the same. Indeed, the wiggles in both 
pressure and acceleration are much smaller in 2D than in ID, for DISPH and SPSPH. In the case of 
SSPH, the wiggle of the pressure is not much different for ID and 2D calculations. The error in the 
acceleration is somewhat smaller in 2D, but this effect is not so drastic as in the case of DISPH and 
SPSPH. Thus, we can expect that the improvement of SPSPH over DISPH is quite large, for realistic, 
multi-dimensional calculations. 


ldim 2dim 




Fig. 13. Results of the one (left) and two (right) dimensional hydrostatic equilibrium test at t = 0. 


5. Discussion 

5.1. Limiting cases for the diffusion constant 

As we discussed in section 2.2, if we set D — 0, our SPSPH is reduced to either SSPH or the 
scheme proposed by Ott & Schnetter (2003). Here, we consider the other limit of D = oo. This means 
we would effectively solve an elliptic equation, instead of the parabolic diffusion equation. Therefore, 
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the pseudo density y would take the same value everywhere. If we set that value to unity, we have 

J2Z b W ab (h a ) = 1, (104) 

b 

from Eq.(46). We can see that, in Eq.(104), Z b is determined purely from the positions of particles. 
Thus, one might think it would give an even better way to formulate SPH equations. We have per¬ 
formed some experiments with this form, so far with little success. It is hard to determine the value 
of Z b in Eq.(104). 

With a matrix Q ab = W ab and vectors R„ = Z a . s a = 1, we can rewrite Eq.(104) into, 

QR = s. (105) 

The condition for the existence of a unique solution is that the matrix Q is regular. In other words, 
all eigenvalues of Q should be nonzero. We introduce two assumptions for simplicity. First, we 
consider the one dimensional case. Second, we assume a uniform particle separation Ax. With these 
assumptions, the eigenvectors q of Q are given by 

q a = e ikx f (106) 

where k is a wavenumber and satisfies 0 < k < ‘In f Ax. Hence, the eigenvalues are 

Yw ab d Kxh - Xa) = T,^Z b W ab e ik ^\ 

b b Z b 

= (107) 

Here, we used the fact that for a uniform particle separation the value of Z b is constant anywhere. 
Eq.(107) means that the fourier transform of W must be nonzero in the range 0 < /;: < 2n/Ax in 
order to have an inverse matrix of Q. However Dehnen & Aly (2012) showed that a fourier transform 
of some kernel functions become zero at 3 k with the uniform particle separation. Moreover, we 
are not sure that those of other kernels such as Wendland functions do not become zero at 3 k for 
inhomogeneous particle distributions. Therefore SPH with functions Eq.(104) is less than successful. 

5.2. Free surface 

The SSPH scheme cannot handle the free surface well, since the density of particles near the 
surface is grossly underestimated. Here we discuss the possibility to extend our method to handle 
free surface. As a simple example, we consider the water surface (like that of sea surface). From the 
physical point of view, the surface of water is not free, but covered by the atmosphere. In other words, 
it is simply the contact discontinuity of water and air. Thus, if we express air as well as water, by SPH 
particles, the surface of the water will be handled properly. This is however impossible with SSPH, 
since it cannot handle large density jumps. However, with our SPSPH, density jump would not cause 
problems. Thus, one solution to the treatment of the free surface is to introduce SPH particles that 
represent thin air. This scheme would work fine for engineering problems in which there actually is 
air. It might also works for problems in planetary science, like the giant-impact simulations. 
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6. Summary 


The SSPH scheme cannot handle discontinuities in density. The reason is that SSPH requires 
that density is positive and spatial continuity even at the discontinuities of density. To solve this 
problem we introduce a new quantity pseudo-density and require the continuity of pseudo-density 
instead of that of the density. Pseudo-density evolve with artificial diffusion for guarantee the positive 
and spatial differentiable. SPSPH can handle the contact discontinuities quite well and has possibility 
for handling the free surface with particles that represent thin air. 
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